home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / kcl / akcl / akcl1615.lha / mp / mpi.c < prev    next >
C/C++ Source or Header  |  1992-02-07  |  15KB  |  645 lines
  1.  
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  4. /*~                                    ~*/
  5. /*~               OPERATIONS DE BASE (NOYAU)            ~*/
  6. /*~                                    ~*/
  7. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  8. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  9. /*   This file was modified by W. Schelter to be suitable for optimization
  10.      and inlining of assembler for maximum speed
  11.  */
  12.  
  13. #include "config.h"
  14. #include "genpari.h"
  15. #include "arith.h"
  16.  
  17. GEN   mulsi(x,y)
  18.      long x;
  19.      GEN y;
  20. { TEMPVARS
  21.   long s=signe(y),ly=lgef(y),i;
  22.   GEN z,zp,yp; ulong hiremainder;
  23.   
  24.   if((!x)||(!s)) return gzero;
  25.   if(x<0) {s= -s;
  26.        x= -x;
  27.        if (x < 0) /* -2^31 */
  28.          {return mulii(stoi(1<<31),y);}
  29.      }
  30.   z=cgeti(ly+1);
  31.   hiremainder=0;
  32.   MP_START_LOW(yp,y,ly);    MP_START_LOW(zp,z,ly+1);
  33.   i = MP_COUNT_LG(ly);
  34.   WHILE_COUNT(--i) { MP_NEXT_UP(zp) = addmul(x,MP_NEXT_UP(yp));}
  35.   if(hiremainder) {MP_NEXT_UP(zp)=hiremainder;
  36.            setlgef(z,ly+1);}
  37.   else {avma+=4;z[1]=z[0]-1;z++;setlgef(z,ly);}
  38.   setsigne(z,s);return z;
  39. }
  40.  
  41. int expi(x)
  42.      GEN x;
  43. {
  44.   long lx=x[1]&0xffff;
  45.   
  46.   return lx==2 ? -8388608 : ((lx-2)<<5)-bfffo(x[2])-1;
  47. }
  48.  
  49. GEN addsi(x,y)
  50.      long x;
  51.      GEN y;
  52. {
  53.   long sx,sy,ly,p,i;
  54.   ulong overflow;  
  55.   GEN z; TEMPVARS
  56.  
  57.   
  58.   if(!x) return icopy(y);
  59.   sy=signe(y);if(!sy) return stoi(x);
  60.   if(x<0) {sx= -1;
  61.        x= -x;
  62.        if (x < 0)  /* x=-2^31 */
  63.          return addii(MOST_NEGS,y);
  64.      } else sx=1;
  65.   ly=lgef(y);
  66.   if(sx==sy)
  67.     {
  68.       p=addll(x,y[ly-1]);
  69.       if(overflow)
  70.     {
  71.       z=cgeti(ly+1);z[ly]=p;
  72.       for(i=ly-1;(i>2)&&(y[i-1]==0xffffffff);i--) z[i]=0;
  73.       if(i>2)
  74.         {
  75.           z[i]=y[i-1]+1;i--;while(i>=3) {z[i]=y[i-1];i--;}
  76.           z[2]=z[1]=z[0]-1;z++;avma+=4;
  77.         }
  78.       else {z[2]=1;z[1]=z[0];}
  79.     }
  80.       else
  81.     {
  82.       z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
  83.     }
  84.       setsigne(z,sx);
  85.     }
  86.   else
  87.     {
  88.       if(ly==3)
  89.     {
  90.       if((ulong)y[2]>(ulong)x)
  91.         {
  92.           z=cgeti(3);z[1]=(sy<<24)+3;z[2]=y[2]-x;return z;
  93.         }
  94.       if(y[2]==x) return gzero;
  95.       z=cgeti(3);z[1]=((-sy)<<24)+3;z[2]=x-y[2];return z;
  96.     }
  97.       p=subll(y[ly-1],x);
  98.       if(overflow)
  99.     {
  100.       z=cgeti(ly);z[ly-1]=p;
  101.       for(i=ly-2;!(y[i]);i--) z[i]=0xffffffff;
  102.       z[i]=y[i]-1;
  103.       if((i>2)||z[i]) {i--;for(;i>=1;i--) z[i]=y[i];}
  104.       else
  105.         {
  106.           z[2]=z[1]=z[0]-1;z++;avma+=4;setsigne(z,sy);
  107.         }
  108.     }
  109.       else
  110.     {
  111.       z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
  112.     }    
  113.     }
  114.   return z;
  115. }
  116.  
  117. GEN addii(x,y)
  118.      GEN x,y;
  119. {
  120.   long sx,sy,sz,lx,ly,i,j,p;
  121.   GEN z,xp,yp,zp,xpp,xhigh; TEMPVARS
  122.   ulong overflow;  
  123.   lx = lgef(x);
  124.   ly = lgef(y);
  125.   if (lx< ly) {z=x;x=y;y=z; sx=lx; lx=ly; ly=sx;} 
  126.        /* lx>=ly so sx==0 ==> sy==0 */
  127.   
  128.   if (0 == (sy=signe(y))) return icopy(x);
  129.   sx = signe(x);
  130.  
  131.   if(sx==sy)
  132.     {
  133.       z=cgeti(lx+1);overflow=0;
  134.       MP_START_LOW(zp,z,lx+1);
  135.       MP_START_LOW(xp,x,lx);
  136.       MP_START_LOW(yp,y,ly);
  137.  
  138. #ifdef QUICK_LOOP
  139.        i = ly - 2;
  140.       QUICK_LOOP(i,ADDXCC);
  141. #else      
  142.       i = MP_COUNT_LG(ly);
  143.       WHILE_COUNT(--i)
  144.     {ADDLLX(MP_NEXT_UP(xp),MP_NEXT_UP(yp),MP_NEXT_UP(zp));}
  145. #endif
  146.  
  147.      
  148.       if(overflow)
  149.     {  GEN xhigh = &MP_HIGH(x,lx);
  150.      again:       
  151.       { GEN xpp = &MP_NEXT_UP(xp);
  152.        if (xpp >= xhigh)
  153.          { if (*xpp == 0xffffffff)
  154.             { MP_NEXT_UP(zp)=0;
  155.               goto again;}
  156.           
  157.          else
  158.            { MP_NEXT_UP(zp) = *xpp + 1;
  159.          while ((xpp = &MP_NEXT_UP(xp)) >= xhigh)
  160.            { MP_NEXT_UP(zp) = *xpp ;}
  161.          z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;}}
  162.        else {z[2]=1;z[1]=x[1]+1;}
  163.     }}
  164.       else
  165.     { j = COUNT(lx - ly);
  166.         WHILE_COUNT( --j)
  167.           { MP_NEXT_UP(zp) = MP_NEXT_UP(xp);}
  168.       z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;
  169.     }
  170.     }
  171.   else
  172.     {
  173.       if(lx==ly)
  174.     /* we have to compare x and y */
  175.     { j = MP_COUNT_LG(lx);
  176.       MP_START_HIGH(xp,x,lx);
  177.       MP_START_HIGH(yp,y,lx);
  178.       WHILE_COUNT(--j)
  179.         { ulong tx = MP_NEXT_DOWN(xp);
  180.           ulong ty = MP_NEXT_DOWN(yp);
  181.           if ( ty > tx)
  182.         {z=x;x=y;y=z;sz=sx;sx=sy;sy=sz;
  183.          goto DIFFER;}
  184.           else
  185.         if ( tx > ty)
  186.           {goto DIFFER;}}
  187.       SAME:     return gzero;  
  188.     DIFFER:;
  189.     }
  190.       
  191.       z=cgeti(lx);overflow=0;
  192.       MP_START_LOW(xp,x,lx);MP_START_LOW(yp,y,ly);MP_START_LOW(zp,z,lx);
  193.       i = MP_COUNT_LG(ly);
  194.  
  195. #ifdef QUICK_LOOP
  196.        i = ly - 2;
  197.       QUICK_LOOP(i,SUBXCC);
  198. #else      
  199.       i = MP_COUNT_LG(ly);
  200.       WHILE_COUNT(--i)
  201.     {SUBLLX(MP_NEXT_UP(xp),MP_NEXT_UP(yp),MP_NEXT_UP(zp));}
  202. #endif
  203.       
  204.       if(overflow)
  205.     { ulong tx ;  
  206.       while((tx=MP_NEXT_UP(xp)) == 0)
  207.         MP_NEXT_UP(zp) = 0xffffffff;
  208.       if (xp >= (xhigh = &MP_HIGH(x,lx)))
  209.         { MP_NEXT_UP(zp) = tx -1;
  210.           while ((xpp = &MP_NEXT_UP(xp)) >= xhigh)
  211.         { MP_NEXT_UP(zp) = *xpp;}}
  212.     }
  213.       else
  214.     { i = COUNT(lx - ly);
  215.       WHILE_COUNT(--i)
  216.         MP_NEXT_UP(zp) = MP_NEXT_UP(xp);
  217.     }
  218.       if(z[2]) z[1]=x[1];
  219.       else
  220.     { zp = &z[3];
  221.       while (*zp ==0){zp++;}  /* x was != y by above */
  222.       zp -= 2;
  223.       i = zp - z;
  224.       zp[1] = (zp[0] = z[0]-i);
  225.       z = zp;
  226.       setsigne(z,sx);
  227.       avma+=(i<<2);
  228.     }
  229.     }
  230.   return z;
  231. }      
  232.  
  233. GEN mulss(x,y)
  234.      long x,y;
  235. {
  236.   long s,p1;
  237.   GEN z; ulong hiremainder;
  238.   
  239.   if((!x)||(!y)) return gzero;
  240.   s=1;
  241.   if(x<0)
  242.     {s= -1;
  243.      x= -x;
  244.      if (x<0)
  245.        return mulsi(y,stoi(x));
  246.    }
  247.   if(y<0) {s= -s;
  248.        y= -y;
  249.        if(y<0)
  250.          return mulsi((s > 0 ? x : -x),ABS_MOST_NEGS);
  251.      }
  252.   p1=mulll(x,y);
  253.   if(hiremainder) {z=cgeti(4);z[2]=hiremainder;z[3]=p1;}
  254.   else {z=cgeti(3);z[2]=p1;}
  255.   z[1]=z[0];setsigne(z,s);return z;
  256. }
  257.  
  258.  
  259. GEN mulii(x,y)
  260.      GEN x,y;
  261. {
  262.   long i,j,lx=lgef(x),ly=lgef(y),sx,sy,lz,p1,p2;
  263.   GEN z; TEMPVARS
  264.     
  265.     GEN zz,yy,zp,xx;
  266.   GEN ylow;
  267.  ulong hiremainder; 
  268.   ulong overflow;
  269.   
  270.   sx=signe(x);if(!sx) return gzero;
  271.   sy=signe(y);if(!sy) return gzero;
  272.   if(sy<0) sx= -sx;
  273.   if(lx>ly) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;}
  274.   lz=lx+ly-2;if(lz>=0x10000) err(muler1);
  275.   z=cgeti(lz);z[1]=z[0];setsigne(z,sx);
  276.  
  277.   MP_START_LOW(xx,x,lx);
  278.   p1 = MP_NEXT_UP(xx);
  279.  
  280.   hiremainder=0;
  281.   i = COUNT(ly-2);
  282.   MP_START_LOW(yy,y,ly);
  283.   MP_START_LOW(zz,z,lz);
  284.  
  285.   WHILE_COUNT (--i)
  286.     { MP_NEXT_UP(zz) = addmul(p1,MP_NEXT_UP(yy));}
  287.  
  288.   MP_NEXT_UP(zz) = hiremainder;
  289.  
  290.   /* restart zz one above bottom */ 
  291.   MP_START_LOW(zz,z,lz); 
  292.   
  293.   MP_START_LOW(ylow,y,ly);
  294.   ly  = COUNT(ly - MP_CODE_WORDS);
  295.   lx -= MP_CODE_WORDS;
  296.   while (--lx > 0)        /* one less iteration first term of x, already used */
  297.     { long tem;
  298.       register long p11;
  299.       p11 = MP_NEXT_UP(xx);
  300.       i = ly;
  301.       yy = ylow;
  302.       zp = &MP_NEXT_UP(zz);    /* *zp = second from low word of z first time through */
  303.       tem = 0;
  304.       /* ZerO is just a 68k kludge to getit to keep 0 in a reg during this loop*/
  305. #undef ZERO
  306. #define ZERO ZerO
  307.       {  int ZerO = 0; 
  308.       WHILE_COUNT(--i)
  309.     { p2 = MP_NEXT_UP(yy);
  310.       p2 = mulul(p2,p11,hiremainder);
  311.         MP_NEXT_UP(zp);
  312.       p2 = add_carry(p2,*zp,hiremainder);
  313.       p2 = add_carry(p2,tem,hiremainder);
  314.       *zp = p2;
  315.       tem = hiremainder;
  316.     }
  317.        }
  318.       MP_NEXT_UP(zp) = hiremainder;
  319.  
  320. #undef ZERO
  321. #define ZERO 0
  322.  
  323.     }
  324.   if(!MP_HIGH(z,lz))
  325.     {                /* shift header one along  decreasing lg and lgef */
  326.       z[2]=z[1]-1;z[1]=z[0]-1;z++;avma+=4;
  327.     }
  328.   return z;
  329. }
  330.  
  331. GEN confrac(x)
  332.      GEN x;
  333. {
  334.   long lx=lg(x),ex= -expo(x)-1,ex1,av=avma,ly,ey;
  335.   long lr,nbdec,k,i,j; ulong hiremainder;
  336.   GEN y,res; TEMPVARS
  337.   
  338.   ey=((lx-2)<<5)+ex;ly=(ey+63)>>5;y=cgeti(ly);ex1=ex>>5; /* 95 dans mp.s faux? */
  339.   for(i=0;i<ex1;i++) y[i]=0;
  340.   ex&=31;
  341.   if(!ex) for(j=2;j<lx;j++) y[i++]=x[j];
  342.   else
  343.     {
  344.       k=0;
  345.       for(j=2;j<lx;j++) {y[i++]=shiftlr(x[j],ex)+k;k=hiremainder;}
  346.       y[ly-2]=k;
  347.     }
  348.   y[ly-1]=0;
  349.   nbdec=ey*0.30103+1;lr=(nbdec+17)/9;res=cgeti(lr);
  350.   *res=nbdec;
  351.   for(j=1;j<lr;j++)
  352.     {
  353.       hiremainder=0;
  354.       for(i=ly-1;i>=0;i--) y[i]=addmul(y[i],1000000000);
  355.       res[j]=hiremainder;
  356.     }
  357.   avma=av;return res;
  358. }
  359.  
  360. /* x/y : uses hiremainder for return */
  361. GEN divss(x,y)
  362.      long x,y;
  363. {
  364.   long p1; 
  365.   
  366.   if(!y) err(diver1);
  367.   if (x == (1<<31))
  368.     return divis(stoi(x),y);
  369.   hiremainder=0;p1=divll((ulong)abs(x),(ulong)abs(y));
  370.   if(y<0) {hiremainder= -((long)hiremainder);p1= -p1;}
  371.   if(x<0) p1= -p1;
  372.   return stoi(p1);
  373. }
  374.  
  375. GEN modss(x,y)
  376.      long x,y;
  377. {
  378.   long y1; ulong hiremainder;
  379.   
  380.   if(!y) err(moder1);
  381.   if (x == (1<<31))
  382.     return modis(stoi(x),y);
  383.   hiremainder=0;divll(abs(x),y1=abs(y));
  384.   if(!hiremainder) return gzero;
  385.   return (((long)hiremainder)<0) ? stoi(y1-hiremainder) : stoi(hiremainder);
  386. }
  387.  
  388. GEN resss(x,y)
  389.      long x,y;
  390. { ulong hiremainder;
  391.   if(!y) err(reser1);
  392.   hiremainder=0;divll(abs(x),abs(y));
  393.   return (y<0) ? stoi(-((long)hiremainder)) : stoi(hiremainder);
  394. }
  395.  
  396. /* uses hiremainder for return */
  397. GEN divsi(x,y)
  398.      long x;
  399.      GEN y;
  400. { 
  401.   long s=signe(y),ly=lgef(y),p1;
  402.  
  403.   if(!s) err(diver2);
  404.   if((!x)||(ly>3)||(y[2]<0)) {hiremainder=x;return gzero;}
  405.   if (x== 1<<31)
  406.     return divii(stoi(x),y);
  407.   hiremainder=0;p1=divll(abs(x),y[2]);
  408.   if(signe(y)<0) {hiremainder= -((long)hiremainder);p1= -p1;}
  409.   if(x<0) p1= -p1;
  410.   return stoi(p1);
  411. }
  412.  
  413. /* this uses the GLOBAL hiremainder to return its remainder
  414.    We cannot make it a local.
  415.  */
  416. GEN divis(y,x)
  417.      long x;
  418.      GEN y;
  419. { ulong hi;
  420.   long s=signe(y),ly=lgef(y),i,d;
  421.   GEN z;
  422.   
  423.   if(!x) err(diver4);
  424.   if(!s) {hiremainder=0;return gzero;}
  425.   if(x<0) {s= -s;x= -x;
  426.        if (x < 0)
  427.          return divii(y,stoi(x));
  428.      } 
  429.   if((ulong)x>(ulong)y[2])
  430.     {
  431.       if(ly==3) {hiremainder=itos(y);return gzero;}
  432.       else {z=cgeti(ly-1);d=1;hi=y[2];}
  433.     }
  434.   else {z=cgeti(ly);d=0;hi=0;}
  435.   for(i=d+2;i<ly;i++) z[i-d]=divul(y[i],x,hi);
  436.   z[1]=z[0];setsigne(z,s);
  437.   hiremainder= (s < 0 ? -((long)hi) : hi) ;
  438.   return z;
  439. }
  440.  
  441. /*
  442. #       entree : x pointe sur i2 de type I (dividende)              
  443. #                y pointe sur i1 de type I (diviseur)               
  444. #                z  contient un pointeur sur le reste si l'on       
  445. #                veut a la fois q et r, 0 si l'on ne veut que le    
  446. #                quotient, -1 si l'on ne veut que le reste          
  447. #       sortie : resultat est  q si celui-ci est attendu, et sinon r 
  448. #                r. z  pointe sur r si q et r sont attendus
  449. #                (toutes les zones sont creees)                     
  450. #       remarque : il s'agit de la 'fausse division' ; le reste est 
  451. #                 du signe du dividende                             
  452. */
  453.  
  454. GEN dvmdii(x,y,z)
  455.      GEN x,y,*z;
  456. {
  457.   long av,lx,ly,lz,i,j,dec,sh,k,k1,sx=signe(x),sy=signe(y);
  458.   long saux,k3,k4,av1,flk4;
  459.   ulong si,qp;
  460.   GEN p1,p2,p3,p4,xp,p1p,p2p,pp; TEMPVARS
  461.   ulong hiremainder,overflow;
  462.   
  463.   if(!sy) err(dvmer1);
  464.   if(!sx)
  465.     {
  466.       if(((long)z==0xffffffff)||((long)z==0)) return gzero;
  467.       *z=gzero;return gzero;
  468.     }
  469.   lx=lgef(x);ly=lgef(y);lz=lx-ly;
  470.   if(lz<0)
  471.     {
  472.       if((long)z==0xffffffff) return icopy(x);
  473.       if(z==0) return gzero;
  474.       *z=icopy(x);return gzero;
  475.     }
  476.   av=avma;if(sx<0) sy= -sy;
  477.   if(ly==3)
  478.     { int lgp1;
  479.       si=y[2];
  480.       MP_START_HIGH(xp,x,lx);
  481.       if (si > (ulong)MP_HIGH(x,lx))
  482.     { lgp1=lx-1; hiremainder= MP_NEXT_DOWN(xp);}
  483.       else { lgp1=lx; hiremainder=0;}
  484.       p1 = cgeti(lgp1); i = MP_COUNT_LG(lgp1);
  485.       MP_START_HIGH(p1p,p1,lgp1);
  486.       WHILE_COUNT(--i) { MP_NEXT_DOWN(p1p) = divll(MP_NEXT_DOWN(xp),si);}
  487.  
  488.       if((long)z==0xffffffff)
  489.     {
  490.       avma=av;if(!hiremainder) return gzero;
  491.       p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;return p2;
  492.     }
  493.       if(lgp1!= 2) {p1[1]=p1[0];setsigne(p1,sy);} else {avma=av;p1=gzero;}
  494.       if(z==0) return p1;
  495.       if(!hiremainder) *z=gzero;
  496.       else {p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;*z=p2;}
  497.       return p1;
  498.     }
  499.   else
  500.     {
  501.       p1=cgeti(lx);
  502.       sh=bfffo(y[2]);
  503.       if(sh)
  504.     { GEN p2p,yp;
  505.       MP_START_HIGH(yp,y,ly);
  506.       p2=cgeti(ly);
  507.       k=shiftl(MP_NEXT_DOWN(yp),sh);
  508.       MP_START_HIGH(p2p,p2,ly);
  509.       i = MP_COUNT_LG(ly-1);
  510.       WHILE_COUNT(--i)
  511.         {
  512.           k1=shiftl(MP_NEXT_DOWN(yp),sh);
  513.           MP_NEXT_DOWN(p2p) = k + hiremainder;
  514.           k = k1;
  515.         }
  516.       MP_NEXT_DOWN(p2p) = k ; k=0;
  517.  
  518.       MP_START_HIGH(xp,x,lx);
  519.       MP_START_HIGH(p1p,p1,lx);
  520.       MP_NEXT_UP(p1p) ;    /* yes go out of range !! */
  521.       i = MP_COUNT_LG(lx);
  522.       WHILE_COUNT (--i)
  523.         { k1 = shiftl(MP_NEXT_DOWN(xp),sh);
  524.           MP_NEXT_DOWN(p1p) = k + hiremainder; k = k1;
  525.         }
  526.       MP_NEXT_DOWN(p1p) = k;
  527.     }
  528.       else {
  529.     MP_START_HIGH(xp,x,lx);
  530.     MP_START_HIGH(p1p,p1,lx);
  531.     MP_NEXT_UP(p1p) ;    /* yes go out of range !! */
  532.     MP_NEXT_DOWN(p1p) = 0;
  533.     j = MP_COUNT_LG(lx);
  534.     WHILE_COUNT (-- j)
  535.       { MP_NEXT_DOWN(p1p) = MP_NEXT_DOWN(xp);}
  536.     p2 = y;}
  537.       si=p2[2];saux=p2[3];
  538.       MP_START_HIGH(p1p,p1,lx); MP_NEXT_UP(p1p) ; /* out of bound */
  539.       i = COUNT(lz+1);
  540.       WHILE_COUNT(--i)       
  541.  
  542.     { GEN pp;
  543.       if(MP_NEXT_DOWN(p1p)==si) 
  544.         {  /* Using fact that next_down does post increment */
  545.         qp=0xffffffff;k=addll(si,*p1p);    
  546.                         
  547.         }
  548.       else
  549.         {
  550.           hiremainder=p1p[-1];qp=divll(*p1p,si);
  551.           overflow=0;k=hiremainder;
  552.         }
  553.       if(!overflow)
  554.         {
  555.           k1=mulll(qp,saux);k3=subll(k1,p1p[1]);k+=overflow;
  556.           flk4=((ulong)hiremainder>(ulong)k);k4=subll(hiremainder,k);
  557.           while(flk4) {qp--;k3=subll(k3,saux);
  558.                k4-=overflow;flk4=((ulong)k4>(ulong)si);
  559.                k4=subll(k4,si);}
  560.         }
  561.       hiremainder=0;
  562.  
  563.       j = MP_COUNT_LG(ly);
  564.       MP_START_LOW(pp,p1p,ly-2);
  565.       MP_START_LOW(p2p,p2,ly);
  566.       WHILE_COUNT(--j)       
  567.         { GEN ppp;
  568.           k1=addmul(qp,MP_NEXT_UP(p2p));
  569.           ppp = &MP_NEXT_UP(pp);
  570.           *ppp =subll(*ppp,k1);hiremainder+=overflow;
  571.         }
  572.       if((ulong)p1p[-1]<(ulong)hiremainder)
  573.         {
  574.           overflow=0;qp--;
  575.           j = MP_COUNT_LG(ly);
  576.           MP_START_LOW(pp,p1p,ly-2);
  577.           MP_START_LOW(p2p,p2,ly);
  578.           WHILE_COUNT(--j){ GEN ppp = &MP_NEXT_UP(pp);
  579.                 ADDLLX(*ppp,MP_NEXT_UP(p2p),*ppp);}
  580.  
  581.         }
  582.       p1p[-1] = qp;
  583.     }
  584.       av1=avma;
  585.       if((long)z!=0xffffffff)
  586.     {ulong lgp3 = lz + 2;
  587.      MP_START_LOW(p1p,p1,lgp3);
  588.      if (p1[1]) {lgp3++;}
  589.        else if (lz==0) sy=0;
  590.      p3 = cgeti(lgp3);
  591.      MP_START_LOW(pp,p3,lgp3);
  592.      j = MP_COUNT_LG(lgp3);
  593.      WHILE_COUNT(--j)
  594.        {MP_NEXT_UP(pp) = MP_NEXT_UP(p1p) ;}
  595.      if(lgp3<3) {p3[1]=2;} else {p3[1]=p3[0];setsigne(p3,sy);}
  596.     }
  597.       if(z!=0)
  598.     {
  599.       for(j=lz+2;(j<lx)&&(!p1[j]);j++);
  600.       if(j==lx) p4=icopy(gzero);
  601.       else
  602.         {
  603.           p4=cgeti(lx-j+2);p4[1]=p4[0];
  604.           if(!sh) for(i=2;j<lx;j++,i++) p4[i]=p1[j];
  605.           else
  606.         {
  607.           hiremainder=0;k1=shiftlr(p1[j++],sh);k=hiremainder;
  608.           if(k1) {p4[2]=k1;dec=1;} else
  609.             {p4[1]=p4[0]-1;p4++;avma+=4;p4[1]=p4[0];dec=0;}
  610.           for(i=2+dec;j<lx;j++,i++)
  611.             {
  612.               p4[i]=shiftlr(p1[j],sh)+k;k=hiremainder;
  613.             }
  614.         }
  615.           setsigne(p4,sx);
  616.         }
  617.     }
  618.       if((long)z==0xffffffff) return gerepile(av,av1,p4);
  619.       if((long)z==0) return gerepile(av,av1,p3);
  620.       dec=lpile(av,av1,0)>>2;*z=p4+dec;return p3+dec;      
  621.     }
  622. }
  623. /* machines which provide an inline version of mulul need
  624.  to provide a function for calls where that inlining can't take place
  625.  */
  626. #ifdef NEED_MULUL3
  627. ulong
  628. mulul3(a,b,c)
  629.    ulong a,b,*c;
  630. { return mulul(a,b,*c);}
  631. #endif
  632.  
  633. #ifdef NEED_DIVUL3
  634. ulong
  635. divul3(a,b,c)
  636.    ulong a,b,*c;
  637. { return divul(a,b,*c);}
  638. #endif
  639.  
  640. /*
  641. ;;- Local variables:
  642. ;;- version-control:t
  643. ;;- End:
  644. */
  645.